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ABSTRACT 

Photons emitted during the epochs of Hydrogen (500 < z < 1600) and Helium recombination 
(1600 < z < 3500 for Hell —>> Hel, 5000 < z < 8000 for Helll —)> Hell) are predicted to 
appear as broad, weak spectral distortions of the Cosmic Microwave Background. We present a 
feasibility study for a ground-based experimental detection of these recombination lines, which 
would provide an observational constraint on the thermal ionization history of the Universe, 
uniquely probing astrophysical cosmology beyond the last scattering surface. We find that an 
octave band in the 2-6 GHz window is optimal for such an experiment, both maximizing signal- 
to-noise ratio and including sufficient line spectral structure. At these frequencies the predicted 
signal appears as an additive quasi-sinusoidal component with amplitude about 8 nK that is 
embedded in a sky spectrum some nine orders of magnitude brighter. We discuss an algorithm to 
detect these tiny spectral fluctuations in the sky spectrum by foreground modeling. We introduce 
a Maximally Smooth function capable of describing the foreground spectrum and distinguishing 
the signal of interest. With Bayesian statistical tests and mock data we estimate that a detection 
of the predicted distortions is possible with 90% confidence by observing for 255 days with an 
array of 128 radiometers using cryogenically cooled state-of-the-art receivers. We conclude that 
detection is in principle feasible in realistic observing times; we propose APSERa—Array of 
Precision Spectrometers for the Epoch of Recombination—a dedicated radio telescope to detect 
these recombination lines. 

Subject headings: Astronomical instrumentation, methods and techniques - Methods: obser¬ 
vational - Cosmic background radiation - Cosmology: observations - Recombination - Radio 
continuum: ISM 


1. Introduction 

The cosmological recombination epoch marks an extended period over which electrons recombine with 
protons and Helium nuclei as cosmological expansion and cooling cause the Universe to transition from a 
fully ionized primordial plasma to a gas of almost completely neutral Hydrogen and Helium atoms. In the 
cosmological concordance model, the cosmological recombination epoch spans redshifts 500 < 2 < 1600 for 
Hydrogen and for Helium recombination the corresponding redshifts are 1600 < £ < 3500 for Hell —> Hel 
and 5000 < z < 8000 for Helll —> Hell. At these redshifts, electrons and photons are tightly coupled through 
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energy exchange via Compton scattering and the thermodynamic temperature of electrons and photons are 
almost exactly equal (at T « 3815{(1 + z)/ 1400} K). During Hydrogen recombination, once even a small 
neutral fraction builds up, photons emitted in any free-bound transition to the ground state would almost 
immediately be re-absorbed by a nearby neutral hydrogen atom, thereby effectively compensating for the 
electron captured (Zel dovich et al.p968 Peebles|[l968 Chlnba fc Sunyaev||2007 ). Hydrogen recombination 
therefore relies on free-bound transitions to excited states followed by a trickle down to the ground state 
with the subsequent emission of numerous recombination photons. During this process, atoms are frequently 
photo-ionized by the huge number of CMB photons and after multiple dissociations and recaptures the 
atoms finally reach the ground state, emitting photons in the Lyman-<a resonance or the 2s-ls two-photon 
continuum (|Zeldovich et al.||1968 Peebles 1968) and, at low frequencies, in transitions among excited states 


(Dubrovich 1975| Rybicki & dell’Antonio 1993 Dnbrovich fc Stolyarov|[l997 ) 


Recombination is expected to be stalled as the ambient radiation temperature in the Lyman-<a transition 
wavelength rises and is balanced in a quasi-static equilibrium with populations in the Is and 2 p states. 
Recombination thus depends on the depletion rate of n = 2 atoms and removal of the excess brightness in 
the ambient Lyman-<a line by Hubble expansion and two photon decay from the 2s state. These processes 
occur at a rate that is ~ 10 7 — 10 8 times slower than the spontaneous transition probability of the Lyman-<a 
and play a vital role in controlling the dynamics of recombination. About 57% of all hydrogen atoms become 
neutral through the 2s-ls two-photon channel ( Chlnba fc Sunyaev||2006a ), which has a vacuum decay rate 


of only A 2 sis — 8.22 s 1 (Goppert-Mayer 


1931 


Breit & Teller 1940; Spitzer & Greenstein 


1951 


Goldman 


1989; Labzowsky et al.|2 005). Consequently, hydrogen recombination is expected to be substantially delayed 
compared to what might be expected assuming equilibrium Saha recombination at the average densities and 
temperatures typical for our Universe. 

As the number of photons per baryon is roughly 1.6 x 10 9 , radiative processes including stimulated 
recombination, induced emission and absorption of photons dominate the populating of atomic levels rather 
than collisions (e.g., Chluba et al. 2007 2010). As the Universe gradually expands and recombination 


proceeds by uncompensated bound-bound and free-bound transitions, the level populations of atomic species 
slowly fall out of equilibrium with the radiation field, which causes departures of the CMB spectrum from 
an ideal Planck form. The spectral lines corresponding to these transitions are predicted to appear as 
redshifted additive deviations to the cosmic microwave background (CMB) spectrum, with most of the 


Hydrogen lines originating from redshift z ~ 1300 — 1400 (Rubino-Martm et al. 2006 Chluba & Sunyaev 


2006a); the observable intensity of these spectral lines is furthermore expected to be uniform on the sky 


and unpolarized. The lines are also substantially broadened owing to the extended recombination time (and 
electron scattering in the case of Hell recombination). At low frequencies, adjacent lines overlap substantially 
and hence the cosmological recombination radiation is expected to manifest itself as spectral ‘ripples’ riding 
on a smooth continuum, which are together an additive component of the extragalactic background light 
extending over radio, microwaves and near IR wavelengths. 

Improvements in our understanding of the recombination history within the framework provided by the 
concordance cosmology, the role of 2s-ls two-photon decay and other factors governing the depletion of the 
n = 2 level population, the importance of the contribution of Helium to the recombination process, progress 
in related atomic physics and, last but not the least important, substantial improvements in computing power 


have all contributed to much improved and more realistic modeling of the epoch in question (see Fendt et al. 
2009; Rub ino-Martm et al.| [201Q, for overview of recombination physics). Detailed calculations as described 


Sunyaev & Chluba (2009) and references therein provide a fairly precise estimate of the spectral ripples 


expected due to the recombination of Hydrogen and Helium. 
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In Fig. [TJ we show the predictions from 0.1 GHz up to 3000 GHz. At high frequencies we can distinguish 


features caused by the Lyman-a line, Balmer-continuum and-<a line, Paschen-a and Brackett-<a lines (Rubino- 


Martin et al. 

2006 Chluba & Sunyaev 2006b; 

Chluba et al. 2007). Helium contributes locally at a typical 

leve} 1 ] of ~ 10% — 30% (Rubino-Martm et al. 

2008 1 ; Sunyaev & Chluba 

2009). Owing to the substantial 


line broadening resulting from the extended period of recombination, at low frequencies the recombination 
radiation appears as (i) a smooth continuum to the radio, microwave and near IR backgrounds, plus (ii) more 
easily distinguishable spectral ripples that has a quasi-periodic frequency dependence. It is the detection of 
these ripply features that we focus on in this work because of its distinctive signature that distinguishes it 
from any other known spectral features in the CMB. 
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Fig. 1.— The additive spectral structure expected in the intensity of the uniform extragalactic background 
light owing to cosmological recombination of Hydrogen and Helium ( Rubino-Martm et al.||2QQ6 Chluba & 


Sunyaev 2006a Chluba et al. 2007 Rubino-Martm et al.| 2008). The wide frequency range is successively 
covered in four panels. The additive spectrum is in units of [mJy sr -1 ] in the panels in the top row and [Jy 
sr -1 ] in the bottom row panels. 


1.1. The importance of an experimental detection of the recombination radiation 

Direct observations of the predicted ripples - the unique fingerprint of the recombination era - would be 
a confirmation of the recombination theory, validating our understanding of the associated atomic physics 
as well as the physical processes in these early times, which in itself is a major motivating factor for experi¬ 
mentally detecting them. A detailed discussion on the dependence of the spectral signatures from the epoch 
of recombination on various contributing factors that affect their strength, shape and position is presented 
in Chluba & Sunyaev (2008a) and Sunyaev & Chluba (2009). We list here some of the major motivations 


for an experimental detection. 


1 Enhancements of features caused by additional feedback effects (Chluba & Sunyaev 


2010) were neglected for now. 




























































-4- 


1. Determining key cosmological parameters. Observing the spectral distortions from the epochs of Hy¬ 
drogen and Helium recombination would, in principle, provide an additional way to determine some 
of the key parameters of the Universe. For instance, the dependence of the predicted recombination 
spectrum on and Tq is shown in Fig. 3 and 4 of Chluba & Sunyaev (2008a). It may be noted 


here that most of the recombination of Hydrogen and Helium and thus the spectral signatures that 
arise from these processes occur over redshifts before the formation of the CMB anisotropies that oc¬ 
curred near the peak of the Thomson visibility function: spectral lines from transitions associated with 
the Hydrogen atoms arise from a redshift range of 1300-1400. Hence detections of the recombination 
lines with increasing accuracy provides a way of constraining the thermal evolution of the universe 
beyond the last scattering surface (see below). This also provides a novel method to measure the pre- 
stellar abundance of Helium and to break parameter degeneracies by combining with CMB anisotropy 
measurements. An illustration of the dependence of the predicted final recombination spectrum on 
contribution from Helll —)> Hell and Hell Hel is given in Fig. 11 in Sunyaev & Chluba (2009). 


2. Probing energy release in the pre-recombination era. Detection of cosmological recombination lines 
allows us to gain a better understanding of the thermal history of our Universe on the basis of the 
different redshifts of Hydrogen and Helium recombination. For instance, it is possible to compute the 
recombination spectrum assuming that the ambient radiation field is a distorted blackbody, where the 
distortion is of y -type ( Chluba fc Sunyaev||2009 ). The //-distortion (Zel dovich fc Sunyaev||1969 ) could 
arise from the blue-shifting of photons by Comptonization in the early universe, due to energy-release 
at z < 50000, when the redistribution of photons in energy by Compton scattering becomes inefficient 
( Bnrigana et al.|[l991 Hu fc SiEk|l993 Chluba fc Sunyaev||2012| ). This distortion is described by the 
Compton-// parameter, y = f aN e cdt. COBE/FIRAS observations place an upper limit of 

y < 1.5 x 10 -5 (95% c.L; |Fixsen et al. 1996; Fixsen 2009). The contributions from Hydrogen and 


Helium to the total recombination spectrum depends both on the value of y parameter and when the 
distortion was created (Chluba & Sun yaev|[2 009). Thus the cosmological recombination radiation may 
allow distinguishing between Compton //-distortions that were caused by energy release before or after 
the epoch of recombination. 


3. Testing recombination physics. Todays most advanced recombination codes (Chlu ba" fc Thomas| [2011 
Ah-Haimoud & Hirata 2011) include many subtle atomic physics and radiative transfer effects (e.g., 


Dnbrovich fc Grachev|2005] Chluba fc Sunyaev 20 06b] [Kholnpenko et al.|2007| Switzer fc Hirata[2 008; 


Chluba & Sunyaev 2008b |Hirata 2008 Grin & Hirata 2010). These calculations ensure that the 
science return from Planck and upcoming CMB experiments is not compromised by inaccuracies in 


the recombination model (Rubino-Martin et al. 2010 Shaw & Chluba 2011). However, the detailed 


dynamics of the recombination process is also reflected in the shape and position of the recombination 
features. Any departures of the recombination spectrum from the theoretical predictions will indicate 


presence of some non-standard process (e.g., Chluba & Sunyaev 2009 Chluba 2010). Thus, by observing 
the recombination radiation we can directly confront our understanding of recombination physics with 
experimental evidence. 


Detecting nK amplitude fluctuations in the extragalactic background brightness, which are indeed a tiny 
perturbation to the orders of magnitude larger sky brightness temperature that arises from the CMB, ex¬ 
tragalactic sources and Galactic emission, is a very challenging problem. However, we opine that receiver 
technology and the ability to produce detector arrays have progressed to the point where it is meaningful 
to look at the practical issues. In this first paper, we study the feasibility of such a detection by modeling 
the sky spectrum as recorded by an ideal instrument and examine whether it is at all possible to recover 







































































- 5 - 


the weak signal embedded in the substantially brighter Galactic and extragalactic foregrounds. We discuss 
methods for fitting to data for the recovery or detection of the faint recombination line spectrum when ob¬ 
served embedded in the substantially brighter cosmic radio background. We also compute optimal observing 
frequencies for the detection and associated signal-to-noise ratio for detection with realistic receivers and a 
purpose-built array of spectral radiometers with due consideration to contribution from atmospheric emission 
further guided by radio frequency interference (RFI) over the band. Some of the challenges are very similar 
to those for the detection of the global 21cm signal (see |Patra et al.|20l 3); however, the recombination ripples 
benefit from their unique frequency dependence, which is hard to mimic by other sources and instrumental 
effects. 


2. The signal and noise level for a detection of the recombination radiation 

In this section, we discuss the frequency dependence of the spectral signal and that of the background 
brightness and additive instrument noise to arrive at estimates for signal-to-noise ratio versus frequency for 
detection of the line structure. The successive peaks of the recombination spectral features correspond to 
the spectral lines arising from bound-bound n (n — 1) electron transitions between adjacent principal 
quantum states of Hydrogen, visible at frequency v ~ 4.7 GHz [n/10] -3 (for emission redshift 2 : ~ 1400). 
If we subtract a low order baseline component from the recombination spectrum in the 1.5-7.0 GHz band, 
the resulting spectral structure is as shown in Fig. [2] As a method of removing the foreground and receiver 
response in an observation with an ideal instrument, if a smooth baseline were to be subtracted from a 
recorded sky spectrum, Fig. [2] is what would be expected as a residual. 



Fig. 2.— The ripply recombination line signal expected to be detected in the 1.5-7.0 GHz band. 

We now address the problem of identifying the most suitable frequency range for the experimental 
detection of recombination lines from the Epoch of Recombination by first considering the relative amplitude 
of the spectral ripples arising from recombination to the noise in the detection. The total system noise is the 
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Fig. 3.— The expected intensity of the recombination signal is shown in black. Our model for the intensity 
equivalent to the minimal total system noise is shown as a red dotted line where the receiver noise is 
assumed to have ideal quantum noise limited performance. Also shown as a blue dashed line is the intensity 
corresponding to system noise when observing with a state of the art cryogenically cooled receiver that has a 
noise temperature of 1 K and in green is the sky intensity when observing with an uncooled receiver assuming 
a noise temperature of 14 K. All intensities are in units of Jy sr -1 (surface brightness units) with the x-axis 
in log(z//GHz). 


additive sum of the spatially-varying sky background, consisting of contributions from discrete and diffuse 
Galactic and extragalactic sources in the sky, the cosmic microwave background and the cosmic infra-red 
background, as well as effects of atmospheric emission, plus receiver noise that is at lowest quantum noise 
associated with zero-point fluctuations. This sets a fundamental limit on the achievable detection sensitivity. 
Later in this section we also consider detection using receivers that have more realistic noise temperatures 
consistent with current technology. We also discuss constraints arising from the requirement that within 
the observing band we need to cover at least a distinctive segment of the quasi-periodic cosmological signal, 
which might provide a unique template specific to cosmological recombination that would be difficult to 
mimic by other astrophysical sources, atmospheric emission, RFI or instrumental effects. 


As a representation of the telescope system noise contribution from additive sky radiation, we construct 
a model spectrum of the background sky brightness all the way from 50 MHz up to 4 THz, extending to just 
beyond the redshifted Lyman-a line that is essentially the highest frequency at which we may expect spectral 
features arising from cosmological recombination. We adopt the model presented in | Subrahmanyan fe; Cowsik 


(2013) for the Galactic and extragalactic contributions to sky brightness and derive sky temperatures towards 
the Galactic pole at 150, 408 and 1420 MHz. A fit of a power-law form to these brightness estimates yields 
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a temperature spectral index of —2.5 and a normalization corresponding to sky brightness of 438 K at 
100 MHz; this model spectrum towards the Galactic pole is adopted to represent the contribution to system 
noise from Galactic emission plus extragalactic discrete sources. We also include a component that is an 


estimate of the far and mid-infrared background; this model is derived from data in Table 47 of Leinert et al. 


(1998). For the receiver noise, we compute the system temperature for two cases: (i) assuming a state of the 
art cryogenically cooled receiver with noise temperature of 1 K above quantum noise (S chleeh et al.| |2012) 
and (ii) that for an uncooled receiver with a more realistic noise temperature of 14 K plus quantum noise 
(Belostotski & Haslett 2007 Witvers et al. 2010). We use the am atmospheric model (Paine 2004) with 
typical conditions at the Chajnantoi^] plateau site to derive the atmospheric opacity versus frequency. The 
system temperature is finally translated to above atmosphere using this opacity so that it may be combined 
with the signal intensity derived above for calculating the signal-to-noise ratio. Our adopted model spectrum 
for the above-atmosphere system temperature is shown in Fig. [3] along with the spectrum of the expected 
recombination signal. 



Fig. 4.— Ratio of spacing between peaks of a transitions from n and n + 2 quantum shells of Hydrogen to 
the nominal frequency versus nominal observing frequency. 

As a criterion for the minimum bandwidth required at any observing frequency we choose a spectral 
window that includes at least two adjacent broad recombination spectral lines that arise from n (n — 1) 
(a) transitions. This minimum bandwidth would correspond to the spacing between peaks of a transitions 
from n and n + 2 shells. Such a spectral window would be expected to include a sufficiently high order of 
the variation in the signal so as to give it a distinctive signature. Fig. [4] shows the ratio of this minimum 
bandwidth to the nominal center frequency versus the center frequency. The detection of Recombination 
Epoch spectral lines clearly requires a larger fractional bandwidth at higher frequencies. A ratio less than 
unity implies that an octave bandwidth at the observing frequency would contain more than two adjacent 
recombination lines, satisfying our criterion. The ratio falls below unity for observing frequencies below about 
18 GHz, which implies that detecting the cosmological line spectrum would require bandwidths exceeding an 


2 https://www.cfa.harvard.edu/spaine/am/cookbook/unix/other_examples/Chajnantor.amc 
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octave if the center frequency were to exceed 18 GHz. Receivers with bandwidths wider than an octave are 
susceptible to self-generated radio frequency interference from harmonics of system clocks and local oscillator 
frequencies; therefore, although sensitivity for a detection would undoubtedly improve with wider observing 
bands, for the ultra-sensitive detection experiment being considered herein it would be unwise to attempt 
detection at frequencies above this value. 



Fig. 5.— Variation in the logarithm of amplitude of the ripples representing the signal from recombination 
versus frequency in log (GHz) units. 


As an estimate of the detectable signal from Hydrogen recombination, we require a measure of the 
peak-to-peak spectral variations versus observing frequency. We first fit a spline through the points in the 
recombination line intensity template that are at frequencies corresponding to the peaks of the n —>• (n — 1) 
Hydrogen transitions. We use a redshift that is a mean value, which we determine by comparing the locations 
of the peaks in the predicted recombination line spectrum with the corresponding locations given by the 
Rydberg formula. Subtracting the spline fit gives a residual recombination line spectrum that represents the 
detectable signal in an observation from which a mean spectral baseline has been subtracted. In Fig. [5] we 
plot the strength of the expected recombination line signal defined as half the peak-to-peak amplitude of 
the baseline subtracted ripples expected in the observed band, versus observing frequency. With frequency, 
the expected ripple amplitude increases by orders of magnitude across the frequency range we consider 
here. It may be noted that although in terms of spectral intensity this ripple signal amplitude increases 
with frequency, as shown above the number of spectral ripples within any octave bandwidth decreases with 
increasing frequency. Thus, a frequency regime that maximizes both aspects is sought. (See Sect. 2.2) 


2.1. Signal-to-noise ratio: a matched filter approach to signal detection 

The detection of spectral lines from cosmological recombination involves first subtracting a baseline from 
the observed spectrum to remove the relatively smoother foregrounds and CMB, after which the residual 
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spectral segment may be examined for the expected spectral ripple using a matched filter. We evaluate 
here the signal-to-noise ratio for a detection method where an estimate of the amplitude of any spectral 
ripple that matches a theoretical template is derived from the measurement. The expected signal, following 
baseline subtraction, is similar to a sinusoidal ripple that has a period that varies systematically across the 
spectral segment. The noise or uncertainty in the estimate of the amplitude of the ripple is derived from the 
measurement errors in the measured intensities in the spectral channels, and from the propagation of errors 
from these measurements to the derived estimate of the amplitude of the signal present in the data. The 
measurement errors in channel data depend on the system noise, which originates from the sky brightness 
and receiver noise, and on the channel bandwidth and integration time. 

We consider here the ideal scenario where the measurement data—the bandpass calibrated spectrum of 
the cosmic radio background from which a baseline has been subtracted—contains a spectral ripple that is 
the same as the expected theoretical template without any other residual contamination. In this case where 
the signal amplitude is estimated using a matched filter derived from the template, we obtain a signal-to- 
noise ratio that is limited purely by the noise present in the spectral channels. This case study admittedly 
represents an optimistic estimate of the signal-to-noise ratio attainable. 

The residual spectrum after baseline subtraction is assumed to have N frequency channels, each of 
bandwidth 5b. The sampled residual sinusoidal ripple, which is also the template of the expected ripply 
recombination lines following baseline subtraction, is given by a k and has an amplitude of Aq. From this 
template of the expected signal we define a weighting function w k , which is defined such that the weighted 
sum of ak is Aq. The matched filter is thus the weighted summation of the observed spectrum and the 
matched filtering yields an estimate of the amplitude S of the ripple from recombination: 

N 

s =akWk (i) 

k=0 

The template representing the expectation and the corresponding weighting function may be scaled appro¬ 
priately to search the data for signatures of the ripple from recombination assuming a range of effective 
redshifts for the recombination lines. 


The weighting function is hence similar to the template in that it has the same ripple form but with an 
amplitude wq. If the signal matches the template, the requirement that S = a k w k ought to equal Aq 

leads to the condition that wq = 2/N. 

The uncertainty in the channel data in the residual spectrum a/ c is 


5a k 


T 

V5b • F 


( 2 ) 


where T is the system temperature and t is the total integration time. Propagating this error in the channel 
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data to the error SS in the estimate of the amplitude of the ripple from recombination: 


SS 


N-l 


N 


y^((5a fc ) 2 |w fc | 2 

k =0 


T 

y/5b • t 


N-l 


N 


y^k/ci 2 




TV 


25b • £ 


The last step above follows from the result that the average value of the square of a sinusoidal wave is half 
the amplitude. Substituting the earlier result that wq = 2/N yields: 


5S = T 


N5b • t 


( 3 ) 


Thus the signal-to-noise ratio in the measurement is given by: 


SNR = is 


A 0 I N5b • t _ A 0 [B~A 


( 4 ) 


where B = N5b is the total bandwidth of the observed spectrum. For this estimate it is assumed that noise 
in the different channels is uncorrelated. 


2.2. Choice of Frequency for detecting ripples from cosmological recombination 

In Fig. [6] we show the estimated signal-to-noise ratio versus observing frequency to guide the choice of 
observing frequency. In creating this plot the signal is assumed to be half the peak-to-peak magnitude of the 
ripple, which is the amplitude shown in Fig. |5j The noise is assumed to be the system temperature (from 
Fig. [ 3 ]) scaled by a factor yj2/B, where B is the octave bandwidth about a nominal central frequency. The 
integration time is assumed to be 1 s. The figure shows the signal-to-noise ratio for the case where the system 
noise corresponds to quantum-noise limited receivers (in green); we also show in red the signal-to-noise ratio 
for a cryogenically cooled receiver with 1 K noise temperature and in blue this ratio for the case of uncooled 
14 K receivers. All three traces are assuming an ideal case where the spectral ripple in the measurement 
data matches the predicted template exactly. 

Detection with maximum signal-to-noise ratio suggests observing in the 1-6 GHz band. If we avoid 
the lower end of this band where Galactic HI and terrestrial and satellite down-link related radio frequency 
interference is substantial, the 2-6 GHz band is suggested. Observing in octave bandwidths within this range 
also satisfies the criterion that at least two cycles of ripples of the recombination line spectrum ought to be 
contained in the observed spectral segment (see Fig. |2|. 


3. A modeling of observations of the cosmological recombination spectrum 

The recombination lines from cosmological recombination are detected by a receiver along with fore¬ 
grounds, which are averaged by the telescope beam over its response pattern on the sky. The averaging over a 
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log (Frequency in GHz) 


Fig. 6.— Signal-to-noise ratio versus nominal observing frequency (in log(GHz) scale) for the detection of 
spectral-line signatures from recombination epoch with an octave bandwidth. The green line shows the ratio 
assuming ideal quantum noise limited receivers; the blue line shows the ratio for the more realistic case of 
receivers with 14 K receiver noise. The red line shows this ratio for a receiver which is cryogenically cooled 
to 1 K. The observing time is assumed to be 1 s. 


multitude of sources with different emission spectra, over the sky and along line of sight, results in a detected 
spectrum that would have an unknown form: even if the emissivity of the gas is a power-law form at every 
location, the spectral index does vary across the sky and along line of sight and, therefore, the averaging by 
the telescope beam would result in an observed spectrum that deviates from a single power law. A similar 
effect, related to the superposition of blackbodies inside the beam of an experiment, causes an inevitable 
y -type distortion of the CMB spectrum (Chluba & Sunyaev 2004); the superposition of blackbodies simply is 
not a blackbody anymore ( Zeldovich et al.||1972| Chluba fc Sunyaev||2004 Stebbins||2007 ). It is necessary to 
detect the presence of the recombination line spectrum in the observed spectrum although the line spectrum 
has amplitude that is almost nine orders of magnitude smaller and the additive foreground is a complex 
spectrum of unknown form! By simulating the sky spectrum as would be observed by an ideal system, this 
section addresses the question of whether such a detection is at all possible. 


We choose an octave band from 3 to 6 GHz as the observing band for the simulations presented here. 
This band also meets our criterion of having at least two recombination spectral lines within the band. In 
this section we describe a code we have developed, which simulates an ideal receiver system observing the 
sky over the identified frequency range, and generates the temperature spectrum of the sky as would be 
produced by such an instrument. A calibration method is included in the pipeline. 


We use all-sky maps at 408 MHz (Haslam et al. 1982), 1420 MHz (Reich 1982 Reich & Reich 1986 
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and at 23 GHz (WMAP science data product^]) as input maps to estimate the combined 
galactic and extragalactic brightness contribution to the sky spectrum between 3 and 6 GHz. The sky 
brightness in the Rayleigh-Jeans limit, at every sky pixel and over frequency channels spaced 10 MHz apart 
is derived from these three input temperature maps using a linear interpolation in log-log space. 


Gorski et al.||2005 


We derive a 408-MHz all sky map from the corresponding raw data product available in the LAMBDA^ 
website. The non-destriped 408-MHz map made with a 0?85 beam that is available in the HEALPix ( |G6rsk i 


et al. 2005) R9 nested ordering scheme is first smoothed with a Gaussian beam of FWHM 31 f60 to obtain a 
map with 1° resolution. Following this we degrade the map to an R8 nested ordering representation. This 
final 408-MHz map used in the simulation is a nested R8 HEALPix map in Galactic coordinates with a 
resolution corresponding to beam FWHM 1° and with temperature in mK brightness units. 


The 1420-MHz all sky map was derived from the data available in the OCEANCOLOP^] website. This 
data is in Kelvin units with a beam size of 0?6 and a pixel size 0?25 and is in J2000 celestial coordinates. We 
perform gridding convolution of the data with a Gaussian of appropriate FWHM to yield an image with final 
resolution of 1°. We then transform the map to Galactic coordinates using the celestial HEALPix coordinate 
pixel listing for the R8 nested ordering scheme. Thus the final all-sky map we use at 1420 MHz has a beam 
FWHM 1° and is in Galactic coordinates ordered in the nested R8 HEALPix scheme. 


The deconvolved form of the WMAP 23 GHz all sky map from the LAMBDA website was smoothed to 
1° resolution and formatted in Galactic coordinates in the HEALPix R8 nested ordering scheme. The image 
units were transformed from thermodynamic to brightness temperature. The image intensities correspond 
to a differential measurement; therefore, the map is uncertain in its zero point, which was arbitrarily set. 
We assume a plane-parallel slab model for the galaxy and use the method adopted in Kogut et al. (2011) 
to estimate this uniform component at 23 GHz using all sky images of the absolute brightness of the sky 
at 150, 408 and 1420 MHz. The CMB monopole, dipole and WMAP ILC images were subtracted from all 
these images and the pixels with substantial contamination from discrete sources were blanked using the 
same blanking image that was used in the WMAP analysis. Pixels in each of the images were binned in 
cosecant Galactic latitude (cosecant (b)) over the range 1.0-4.0. The run of median pixel intensity versus 
mean cosecant (b) were examined in log-log coordinates. A straight line fit to this plot yields an estimate 
of the mean extragalactic background at each frequency. It may be noted here that we have restricted this 
analysis to the southern Galactic hemisphere for consistency with the WMAP analysis. For the case of the 
23 GHz image, the fit yields an intercept of OK with an accuracy of 1 part in 10 9 as expected for an image 
has been constructed to have its intercept arbitrarily set to zero. The intercepts estimated for the lower 
frequency images were extrapolated using a polynomial fit to get the ‘missing’ uniform background in the 
23 GHz image. The resulting value of 493 fiK was added as a constant to the 23 GHz image. 

The three maps representing the Galactic and extragalactic radio emission (minus CMB) at 408 MHz, 
1420 MHz and 23 GHz were interpolated at every image pixel separately to estimate sky temperatures at 
frequencies between 3 and 6 GHz. We interpolate using a first order polynomial in log-brightness-temperature 
versus log-frequency space, considering the spectrum at every pixel to be a smooth power-law in linear space. 
The total sky brightness toward any sky pixel and frequency is estimated as the sum of the galactic and 
extragalactic sky brightness derived from this polynomial interpolation, plus a uniform cosmic microwave 


was arbitrarily set 


Kogut et al. (2011 


3 WMAP Science Team 
4 http://lambda.gsfc. nasa.gov/ 

5 http://oceancolor.gsfc.nasa.gov/AQUARIUS/DinnatEtA12010/ 
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Fig. 7.— Sample synthetic spectrum recorded by an ideal receiver system discussed in the text. The receiver 
is assumed to be observing with an antenna with a cos 2 (ZA) beam pattern pointed at the zenith. 


background brightness computed using the Planck formula, plus a uniform component corresponding to the 
weak cosmological recombination line spectrum, as shown in Fig. [l] 


The observing system is assumed to have a single antenna element with a frequency independent radia¬ 
tion pattern that has a cos 2 (ZA) form, where ZA denotes zenith angle. The antenna bore sight (the optical 
axis the antenna) is assumed to be static on the ground at the observing site and directed towards zenith 
at all times. We neglect polarization in the sky intensity as well as the telescope response. Separately, the 
spectral radiometer is assumed to be a correlation spectrometer and hence does not respond to receiver noise 
(see Patra et al. 2013). The spectral response of the observing system at any instant is then simply the 
weighted sum of spectra towards every pixel in the sky where the weighting is by the beam power pattern 
towards that sky direction. As stated above, the spectrum towards any sky pixel is not expected to be a 
simple power law because it is an average along the line of sight of regions that might have different spectral 
indices. The weighted average spectrum is an average over the sky and along the line of sight and hence 
even if each emitting region is of power law form the average observed spectrum may have a complex form 
that is unknown. 


Our simulation code has the geographic location of the observing site and LST as free parameters. Our 
code can be easily modified to change the antenna beam pattern as well. The generation of spectra takes 
into account effects of atmospheric refraction and precession; astronomical aberration is negligible even in its 
severest form for the problem considered here and is dropped from our calculations. Calibration is assumed 
to be done by recording spectra with a hot (373.0 K) and separately a cold (273.0 K) load on the antenna 
and dividing spectra recorded on the sky by the difference between spectra recorded with the hot and cold 
loads. Bandpass calibration of sufficiently high precision is assumed. Our code generates mock calibrated 
spectra over time whose mean temperature varies as the sky and Galactic Plane drift across the telescope 
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Frequency (GHz) 

Fig. 8.— A time-frequency plot showing the synthetic sky spectra as observed over 24 hours. One can see 
the galaxy rise and set as the sky drifts over the instrument with each spectrum recorded one minute apart. 

beam. A sample spectrum of the sky is shown in Fig. [7| Figure [8] shows a waterfall plot to illustrate the 
variation over 24 hours in the recorded spectra observed by this ideal instrument as the sky drifts across the 
antenna. 

The synthetic sky spectra, whose amplitude is of the order of a few K, contains the recombination line 
spectrum as a small additive component and is representative of mock observations made with a correlation 
spectral radiometer that does not respond to receiver noise. The challenge is to distinguish the weak signal 
corresponding to the epoch of cosmological recombination, which is buried in the observed spectrum as a 
broadband quasi-periodic sinusoid with peak-to-peak amplitude of order ~ 10 nK. Not only is the signal 
a tiny fraction of the total sky spectrum, the recombination line spectrum is an additive component of a 
foreground spectrum whose functional form is complex and unknown. In the next section we discuss methods 
to detect the cosmological recombination spectrum and challenges therein. 


4. The detection of signatures of cosmological recombination 

In the previous section we have described the generation of a synthetic spectrum of the radio sky between 
3 and 6 GHz, as would be observed by an ideal instrument. The signal from cosmological recombination 
is expected to appear as quasi-periodic ripples, some 9 orders of magnitude smaller than the galactic and 
extragalactic foreground spectrum in which it is additively concealed. The discernment of the recombination 


Temperature (K) 
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signal in such a total spectrum, even under ideal conditions, is pivotal in answering the question of whether 
an experimental detection is indeed possible. Below we discuss the challenges involved in and describe a 
possible method for the detection of the cosmological recombination signal in synthetic sky spectra. 

The brightness of the sky as measured by a total power radiometer is the sum of the brightness contri¬ 
butions from all the discrete and diffuse Galactic and extragalactic sources that lie within its beam, along 
the line of sight and across the sky, including the cosmic microwave background and other cosmological 
emissions. The precision with which the foreground needs to be modeled so that a subtraction of the model 
might reveal the cosmological recombination spectrum is extreme. Although the discrete sources lying in the 
beam might have well-measured spectra and the diffuse sky radiation might have all-sky maps at multiple 
frequencies, the functional form of the final cumulative spectrum of the foreground is not known a priori to 
the required accuracy. The foreground has necessarily to be fit to the measurement, which implies that the 
functional form used for the fit to the foreground needs to be of a form that accurately fits to the foreground 
without also fitting to the recombination spectral features. Only then may a fitting of a foreground model to 
the measurement set and its subsequent subtraction reveal the recombination line features in the residual. 

The thermal and non-thermal processes (e.g., synchrotron and free-free emission) that contribute to 
the foregrounds have spectra that are smooth because they arise from impulsive emissions in time domain. 
The electron energy spectra in astrophysical objects and diffuse matter are also believed to be smooth over 
decades in energy space. Hence the cumulative spectra in thermal and non-thermal emissions are expected 
to be featureless over octave bandwidths and, therefore, ought to be distinguishable in principle from the 
cosmological recombination radiation. 

The problem of recovering broad and weak spectral deviations that are buried in a total sky spectrum 
that is several order of magnitudes brighter is not unique to cosmological hydrogen recombination lines. 
All-sky (also referred to as global) spectral deviations predicted to arise from the Epoch of Reionization 
(EoR) also face a similar challenge in that broad spectral features of about 10-100 mK brightness are 
buried in foregrounds of several hundred K, the foreground being ^ 5 orders of magnitude brighter.Thus the 
importance of a method to effectively fit a smooth functional form to the combined Galactic and extragalactic 
foreground over 1-2 octaves of bandwidth, in which the signal of interest is itself not lost in the process, 
cannot be understated. 


In the context of detecting global EoR signal in the all-sky radio background spectrum, a polynomial 
functional form of high order has been adopted in the literature as the analytic function to fit to the 
foreground (Bowman & Ro gersp OlO). However, broad spectral lines of interest that are present in observed 
spectra as additive components may also be absorbed in the polynomial fit, particularly if the polynomial 
is of high order, such that the residual would have little trace of the cosmological signal. As an illustration, 
on fitting a sample spectrum generated by our simulation with an eighth-order polynomial in log-frequency 
versus log-temperature space without any constraints on the nature of the polynomial, and subtracting the 
fit spectrum from the original spectrum, we obtain a residual as shown in Fig. [9] While one might hope that 
the polynomial fit to the foreground would leave the recombination line signal untouched, comparison with 
Fig- i shows clearly that this is not the case. The ripples from cosmological recombination, the very signal 
that we are looking to detect, has been eliminated by the fitting process. 


What is needed is a careful choice of the functional form that would leave the cosmological signals of 
interest almost wholly in the residue. Also needed is a method of successive approximation that would model 
the foreground to the required accuracy so that the faint cosmological signal might dominate the residue. 
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Fig. 9.— The residual obtained on subtracting a sample sky spectrum by an eight order polynomial fitting 
function in log-frequency and log-temperature space. While one would expect to see signatures of the spectral 
lines arising from cosmological recombination in the residual at the nK level, the residual is dominated only 
by thermal noise with the lines themselves having been absorbed by the fitting function. 


4.1. A complete monotone approach to foreground modeling 

The CMB component has a well defined functional form. The fitting function /f g (z/) for the foreground 
may be defined to be an analytic function over the frequency range of interest, so that the model does not 
admit discontinuities and its derivatives are defined. The functional form describing the foreground model is 
expected to be ‘smooth’ in the sense that it must not be able to fit to additive signals that have the sinusoidal 
structure expected of the recombination spectrum. 

A potentially useful functional form is that of a completely monotonic function. Mathematically, a 
function f{pc) is said to be a complete monotone if for all values of x in the interval 0<x<oo, (—l) n x 
d ^ > 0 for every integer n > 0. An example of a complete monotone function is f(x) = l/(a + bx ) c , 
where a > 0, b > 0 and c > 0. This implies that power-law form spectra with negative spectral indices 
are completely monotonic functions. It is also known that if f(x) and g{x) are completely monotonic, then 
af(x ) -\-bg(x) where a and b are non-negative constants is also completely monotonic, which implies that the 
sum of power law spectra with negative indices would also be completely monotonic. A function f{x) may 
also be completely monotonic over finite range a < x < b where a > 0 and b > 0: this would be the case if 
the condition (—l) n x d > 0 for every integer n > 0 holds over the range a < x < b. The mathematical 
definition suggests that if /f g (V), where v is frequency, is the foreground model over a certain bandwidth, 
then we may adopt a functional form for /f g (z/) that is a complete monotone within the bandwidth of interest. 

Successive approximation of an analytic functional form to data may be done using a Taylor approxi¬ 
mation. The Taylor series of an analytic function always converges about every point in its domain. We may 
represent /f g (V) as a polynomial whose coefficients are constrained so that /f g (V) is completely monotonic 
and expand the function as a Taylor polynomial by successively estimating the coefficients of the polynomial, 
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stopping when the degree of the polynomial is such that it is a sufficiently good fit to the data (i.e. the 
observed sky spectrum) and the residual is dominated by the embedded cosmological recombination lines. 
It may be noted here that increasing further the order of the polynomial, whose coefficients are constrained 
such that the polynomial is a complete monotone, will only yield vanishingly smaller coefficients for the 
higher order terms. Since the constraint on the polynomial forces it to remain smooth, once the recombi¬ 
nation lines dominate the residual the residual is not minimized further by the introduction of higher order 
terms. Thus we converge to a smooth completely monotonic approximation to the foreground that leaves 
the recombination line spectral structure as a residual to the fit, which no longer changes significantly with 
increasing the order of the completely monotonic polynomial. 

/f g (z/) may be modeled as a completely monotonic (CM) polynomial in brightness temperature versus 
frequency space. More useful is a modeling in temperature versus log-frequency space, where the number of 
terms in the Taylor approximation would be reduced. The Taylor expansion may be conveniently computed 
using the lowest frequency in the range of interest as the reference value. 

A CM polynomial of arbitrary order n, whose constant term ao is left unconstrained so as to allow 
arbitrary vertical translations, may be written in the form 

n n—i 

f(x) = a 0 + ^"(-l)'(.r - x 0 y{^2 a i+ jCj +J (x m - x 0 ) J }, (5) 

i=l j=0 

where C denotes the binomial coefficient n\/{k\(n — &)!}. An example illustrating the algorithm to construct 
such a CM function and thereby arriving at the functional form presented in equation [5] is described in 
Appendix [Aj 

A disadvantage of adopting a CM functional form is that the data is fitted in brightness temperature 
versus log-frequency space, where the foreground is expected to be CM. A functional form that is smooth in 
log-temperature versus log-frequency space is preferred, since a lower order polynomial would be sufficient 
to describe the foreground. Lower numbers of parameters makes for more robust fitting with less likelihood 
that optimization algorithms get trapped in local minima. 


4.2. Modeling the foreground as a Maximally Smooth function: a smooth polynomial that has 

no zero-crossings in derivatives 

We now consider polynomial functional forms to describe the smooth foreground in log-temperature 
versus log-frequency space. The first approximation in this parametrization is a straight line representing 
the mean spectral index of the sky region. Curvature in the mean spectrum may be represented, to lowest 
order, by adopting a parabolic form for the spectrum in this log space. This form has a constant second 
order derivative. 

If we improve upon the modeling of the spectral curvature by adding a cubic term to the polynomial, 
we may constrain the model to be smooth and without inflections by requiring that the second derivative 
has no zero crossings in the domain. If we considering a polynomial model of the form 

f(x) = Po + P!(x - x 0 ) + p 2 (x - x 0 ) 2 + p 3 (x - x 0 ) 3 , (6) 

we would require that second derivative 

= (2!/0!)p 2 + (3!/l!)p3(* - *o) 


(7) 
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has no zero crossings in the domain. 

If we wish to improve, further, the modeling of the spectral curvature we may represent the foreground 
spectrum by a fourth order polynomial of the form: 

f(x) = Po+Pi(x - Xo) +P 2 (x - xo) 2 +Pz(x - xo) 3 +Pa{x - x 0 ) 4 . (8) 

To constrain this polynomial to be smooth and without inflection points embedded, we now require that the 
second derivative 

= (21/0 \) P2 + (3!/1 !)p3 (* - xo) + (4!/2 \)p 4 (x - x 0 ) 2 (9) 

has no zero crossings. Additionally, we also require that the third derivative 

= (3!/0!)p 3 + (4!/l!)p 4 (x - x 0 ) (10) 

has no zero crossings. 

In general, we may model the foreground by an n th -order polynomial and require that all derivatives of 
order 2 and higher have no zero crossings within the domain of interest. This is implemented by computing 
the functions 

d m fix) 

, - = {m\/0\}pm + {(m + l)!/l!}p m +i(x - x 0 ) + {(m + 2)!/2 \}p m+2 (x - x 0 ) 2 + 

ax 171 

{(m + 3)!/3!}p m+3 (a; — xo) 3 + . + {n!/(n - m)\}p n (x - x 0 ) n ~ m (11) 


d m /(x) 

dx m 


n—m 

E 

i =0 


{(m + i)!/i!}p m+i (x - x 0 y 


( 12 ) 


for all m in the range 2, 3, 4, ...., (n — 1) and constraining the polynomial coefficients pj so that there are no 
zero crossings within the domain for any of these functions. We call the polynomial functions that satisfies 
these constraints Maximally Smooth functions. They will not have ripples embedded. 


We model the foreground in log-temperature versus log-frequency space using a Taylor series expansion 
about the lowest frequency logio(^o) i n our band. The polynomial is written in terms of powers of log 10 (z//z/o). 
If y(x) is the polynomial describing the foreground in log-space, then 10 y ^ is added to a term that describes 
the CMB spectrum and another term that models the recombination line spectrum to get a model for the 
total spectrum. 


We model the recombination line spectrum component in the mock data using a single scaling parameter: 
the recombination line component is considered to be this scale factor times a template of the spectral 
ripple that is nominally expected to be present in the observation. First, the scale factor is in itself a 
quantitative measure of the recombination template present in the total spectrum. A scale factor close to 
unity indicates that the predicted template is present in the spectrum where as a small value indicates the 
absence of such a spectral signature. The distribution allowed for this scaling parameter by the goodness of 
fit yields the confidence in the detection of the predicted recombination template. Secondly, by modeling the 
recombination lines using a scaled template allows the residuals of the optimization to approach measurement 
noise if the theoretical framework that led to the template is correct. 

The analytic fitting function in its final form, which includes a Maximally Smooth form for the foreground 
modeling, is given by 

his\ 

t) 


/ (e^o - l) +PiT iec (u) + lo£r=o[i°Sio(^o)]W 3 



(13) 
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where po corresponds to the CMB temperature, pi is the scale factor that multiplies the recombination line 
template T rec (u ) and P 2 through p 2 +n are the coefficients of the terms in the n th -order Maximally Smooth 
polynomial that models the foreground. Fig. [2] shows the predicted recombination line signature that is 
expected to be present in the synthetic spectrum as an additive component; as stated earlier, this template 
has been derived from the predicted recombination line spectrum by subtracting a baseline of low order and 
is what would be expected as a residual if a smooth baseline were to be subtracted from an observation 
as a method of removing the foreground. Segments of this template are what form the recombination line 
template T rec (z/). 


We use the downhill simplex (Nelder & Mead[l965) optimization algorithm to iteratively fit this model 
to the synthetic spectrum adopting a successive approximation strategy. We start by fitting to a function 
that has four coefficients, one for the CMB temperature, one for the amplitude of the recombination line 
spectrum and two that describe a first order model for the foreground. We successively include more terms 
one by one giving an initial guess of zero for each new coefficient. 


Adopting the method discussed in Section 3, we have generated synthetic spectra over three octave 
bands — 2.0-4.0, 2.5-5.0 and 3.0-6.0 GHz — spaced uniformly over 24 hr in LST to test the robustness of 
this modeling algorithm. In this test, the noise in the spectrum was kept small so that the residual spectrum 
might reveal the recombination line structure obviously, if recovered successfully. As expected, the final fit 
parameters vary for different realizations of the sky spectrum corresponding to different observation times. 
Nevertheless, the algorithm did indeed converge for all LSTs and yielded a best fit scaling factor close to 
unity at all times. An 8 th order Maximally Smooth polynomial was used to model the foreground. The 
residuals were consistent with the measurement noise that had been added to the synthetic spectra and the 
foregrounds were indeed successfully modeled as 8 th -order smooth polynomials, without also fitting to the 
embedded recombination line spectrum. 



Fig. 10.— Foreground model fit to a synthetic spectrum that represents a mock observation is shown in 
blue. The residual, scaled by a fraction of the inverse of the chi-square of fit and and added back to the 
foreground model is shown in red. 
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In Fig. [lO] we show a sample fit to a mock observation in the 3-6 GHz band. The model-fit CMB plus 
Maximally Smooth foreground is shown along with a second trace that is constructed by computing the 
residuals to this fit and adding the residuals to the model with a large rescaling. The plot shows that the 
model does indeed not fit to the recombination line ripple, which is present in the synthetic spectrum, and 
does leave the ripple as a residual. The model fitting algorithm based on the Maximally Smooth polynomials 
we have constructed is capable of separation of the foreground from embedded ripple, which is nine orders 
of magnitude smaller. The residual alone is shown in Fig. El which can be compared with Fig. [9] where an 
8 th -order polynomial was used without any constraints to fit to the synthetic spectrum. 



Fig. 11.— The residual on subtracting the foreground plus CMB model-fit from the synthetic spectrum. 





Fig. 12.— Residuals in the 2.0-4.0, 2.5-5.0 and 3-6 GHz bands following fitting and subtracting Maximally 
Smooth foreground models to synthetic sky spectra. The mock observations are over a 24 hour period, with 
spectra spaced 2 hours apart, and represent observations made by an ideal instrument as discussed in the 
text. Shown as a red dotted line is the recombination line template from the theoretical predictions; the 
overlaid colored dashed lines represent the recovered residuals from different realisations of synthetic sky 
spectra. 
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In Fig. [12] we show the residuals to the fit for a smooth foreground plus CMB to synthetic spectra for 
mock observations in different frequency ranges and distributed over LST. The algorithm does recover the 
recombination ripple in all cases, giving confidence in the robustness of the method. 

These residuals represent a detection of the cosmological recombination lines. The recovery demon¬ 
strates that the smooth functional form proposed here does not also fit to the nanoscale embedded ripple 
representing the recombination line signature, despite the functional form being allowed to be of arbitrarily 
high order. The recovery demonstrates that it is indeed possible to model fit a smooth functional form— 
e.g., the Maximally Smooth functional—to observations and model the completely monotonic foreground—of 
unknown functional form—to an accuracy better than 10 nK! 


5. Confidence in Detection 

In the previous section, we presented an algorithm by which we recover cosmological recombination 
lines of amplitude ~10 nK that are a tiny additive part of a sky spectrum which is nine orders of magnitude 
brighter. In this section, we first examine the probability distribution for the detected amplitude of the 
spectral ripple. We then examine the confidence in detection of the expected ripple and the confidence with 
which we may reject false positives, using a Bayes factor approach. These lead to inference of the observing 
time needed for detections with different confidence levels. 


5.1. Detection based on the fit amplitude of the recombination spectral ripple 


An approach to determining whether or not a given sky spectrum contains cosmological Hydrogen 
recombination lines is to jointly fit to the spectrum a model composed of a Planck spectral component for 
the CMB, a Maximally Smooth polynomial accounting for radiation from point and extended sources in 
the foreground, and the recombination line template scaled by a factor, which is a variable for the fitting. 
The functional form for such a model is given by Equation [l3j We expect that the fit value of p \, which 
represents the normalized amplitude of the recombination line ripple, would take on a value of unity if the 
theory is valid, and a value zero if no lines exist. We choose 0.5 as the threshold between a null detection 
and a positive detection, with values smaller than 0.5 assumed to indicate a null detection and greater than 
0.5 deemed to be indicating a positive detection. 


To determine the confidence in detection using such an approach, as well as the probability of a false 
positive, we synthesize two types of mock observations, one in which the theoretically expected cosmolog¬ 
ical recombination lines are added and another without. Forty spectra of each type were generated with 
independent thermal noise. Below we adopt the notation that data and associated terms for spectra with re¬ 
combination lines present in them would be referred to as data set ‘a’ and those without recombination lines 
as L b\ We fit each of the synthetic sky spectra separately with the mathematical model given by Equation p~3] 
using the successive approximation approach, optimizing all the parameters using the Nelder & Meadj(1965) 
algorithm to minimize the Chi-squared difference between the model and mock data, each time increasing 
the degree of the Maximally Smooth polynomial in Equation [13] by unity until the root mean square of the 
residual saturates. 


To sample the distribution in scaling factor pi and thus derive the confidence for detection using a 
threshold of 0.5 for the scaling factor, we adopt a Markov Chain Monte Carlo (MCMC) analysis. Specifi- 
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cally, we adopt the EMCEE package (Foreman-Mackey et al. 2013), which is a Python-scripting-language 
implementation of an affine invariant ensemble sampler as proposed by Goodman & Weare (2010). The 
emcee package wins over traditional implementation of MCMC samplers in its high computational efficiency 
in generating statistically independent samples from the posterior probability distribution function (PDF). 
Another major advantage of emcee is that it requires tuning of only a couple of parameters by hand in an 
N -dimensional space as opposed to N 2 values in other traditional methods. While data driven algorithms 
do exist to arrive at optimum guesses for initial parameters, for conventional methods this comes at the 
computational cost of lengthy burn-in chains. 


Chi-squared minimization of the model T{y) leads to optimum values for the parameters po representing 
the CMB temperature, p\ representing the scaling factor, and the polynomial coefficients P 2 , P 3 , P 4 , P 5 and 
Pq (for a A th order Maximally Smooth polynomial). We then invoke the EMCEE sampler with 30 walkers 
each of which generates a Markov sampling chain for every parameter. The different chains are initialized at 
random locations close to the optimum location in the N = 7 dimensional parameter space. The EMCEE 
ensemble sampler operates on the log posterior probability distribution function that we define to be: 

ln(Probability) = + ln(27r<r 2 )^) , (14) 

where y a ^ is the synthetic sky spectrum with or without recombination lines. We first run a short burn-in of 
1000 steps, discard these values and run a second burn-in of 6000 steps. We ensure that length of the second 
burn-in is at least ten times the correlation length, then discard these as well and run a final useful 500 steps. 
With 30 walkers exploring 500 steps each we generate a distribution of 15000 samples per parameter per 
synthetic sky spectrum. Marginalizing over the nuisance parameters, which in our case are all parameters 
with the exception of the scaling factor, we have for 40 sky spectra a distribution of 600,000 scaling factors 
in all. On running the EMCEE ensemble sampler on data sets ‘a’ and L b’ > as described above we have two 
distributions of scaling factors with 600,000 samples each. 


We generate multiple pairs of data sets ‘a’ and ‘6’ with different noise integration times; as a guide 
we have adopted times corresponding to different signal-to-noise ratios as defined by Equation [4j For each 
integration time, the fraction of samples of the scaling factor from the Markov chain corresponding to data 
set L a > that have values exceeding the threshold of 0.5 gives an estimate of the probability of detection. 
Similarly the fraction of samples corresponding to data set ‘5’ below 0.5 is an estimate of the probability of 
rejecting a false positive. These are depicted in Fig. [l3j 

Detection of the recombination ripple with 68% confidence using uncooled receivers requires observing 
with about 210 x 10 8 antenna secs, or equivalently about 240000 antenna days. An array of 128 singly- 
polarized total-power spectrometer elements with uncooled receivers would require about ^ 5 years for 
detection with this confidence and based on SF threshold. However, using cryogenically cooled state-of- 
the art receivers, the observing time for detection with 68% confidence is ^ 10 x 10 8 antenna secs or, 
equivalently, 92 observing days with an array of 128 spectral radiometers deploying such cooled receivers. A 
90% confidence detection requires observing with 430 days with such an array, and 95% confidence detection 
requires increasing the observing time for such an array to about 625 days. 


5.2. Detection based on Bayes factor (BF) 

Consider some data D and a choice between model Mi characterised by a set of parameters 6 \ and M 2 
characterised by a set of parameters 62 , to describe the data. The Bayes factor is given by the ratio of the 
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Fig. 13.— The confidence with which scaling factor (SF) above threshold 0.5 signifies a detection of cos¬ 
mological recombination lines (shown as a blue solid line with scale on the left of the plot), and percentage 
likelihood of false positives appearing as a detection for this threshold (shown as a red dotted line with scale 
on the right side of the plot). The confidence values and percentage likelihoods of false positives are shown 
versus integration time assuming cryogenic cooled state-of-the-art 1 K receivers and a 128-element array of 
precision total-power spectrometers. 


probability of the data D given the model Mi to the probability of D given M 2 : 

P{D\M X ) _ J P(e 1 \M 1 )P(D\e 1 ,M 1 )d9 1 


BF = 


p{d\m 2 ) f p(e 2 \M 2 )p(D\e2,M 2 )de 2 ' 


(15) 


An alternative expression of the Bayes factor is in terms of posterior and prior odds. Rewriting the proba¬ 
bilities associated with the two models as odds, we obtain: 


Posterior odds = BF x Prior odds, 

PiMAD) _ P(Mi) 
or that : ; = BF x v ' 


P(M 2 \D) 


p{M 2 y 


(16) 


which leads to : 


Bp = P(M 1 |D) x P(M 2 ) 


P(M 2 \D) " P(Mi)' 

When the prior odds are 1:1 the Bayes factor reduces to the likelihood ratio. 


We once again adopt the notation used in Section [5TT] so that symbols representing data and associated 
terms for mock observations with the recombination lines present in them would be given subscript V and 
those without recombination lines V- Terms representing the null hypothesis would have subscript ‘o’ and 
alternative hypothesis V- In our problem of estimating the confidence in the detection of recombination 
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line ripples in observed spectra, we consider two hypotheses. The first hypothesis is that the spectrum does 
not contain recombination lines of any detectable amplitude; z.e., within errors the mock spectrum can be 
adequately modeled by the total of contributions from the cosmic microwave background and a Maximally 
Smooth function representing the cumulative spectra from diffuse and point sources in the foreground. We 
refer to this as the null hypothesis H 0 . The second hypothesis, which we refer to as the alternative hypothesis 
# 2 , is that the observed sky spectrum includes the cosmological Hydrogen recombination lines as an additive 
component along with contributions from the CMB and foregrounds. Hypothesis H 2 is that the spectrum in 
question does contain the recombination lines exactly as given by the theoretical expectations. The Bayes 
factor for a comparison between these alternate hypotheses is given by: 


PlP|g 2 ) P(H 2 ) 
P(D\H 0 ) P(H 0 y 


(17) 


where P(D\H 2 ), P(D\Ho ) are the respective likelihoods of the alternative and null hypotheses and P(i^ 2 ), 
P(Hq) are the prior probabilities of the alternative and null hypotheses. Since in this comparison we have 
no a priori preference or bias towards the presence or absence of the recombination ripple in data, we may 
assume equal prior probabilities and hence the Bayes factor reduces to the likelihood ratio given by: 


BF = 


P(D\H 2 ) 

P(D\H 0 y 


(18) 


We now proceed to compute the likelihood functions for the above two hypotheses. We generate data 
set ‘a’, a set of 100 mock observations with independent thermal noise and with each containing the recom¬ 
bination lines as an additive component. We then fit each of these spectra independently with two different 
models. The first model is the null hypothesis Hq with a blackbody function to model the CMB compo¬ 
nent and a Maximally Smooth polynomial to model the foreground component of the sky spectrum. Hq 
expects that these two components completely describe the spectrum and that the residual on subtracting 
this model-fit from the data must be Gaussian random with zero mean and standard deviation given by the 
measurement noise. We compute the likelihood of obtaining the data given Hq by: 


P(D a \H 0 ) 


N 


n 


_ y vf. .<i(i [d] 

2ag 


(19) 


where N is the number of independent points across the spectrum and y re so[i\ is the residual spectrum 
following subtraction of the model corresponding to the null hypothesis. 

The variance of the measurement noise, Og, is estimated from the data itself. This variance is assumed 
to be half the Allan variance (AY) of the residual, which is given by: 


N-l 

AV = Yi (: Vres[i + 1] - 2/resH) 2 , (20) 

i= 1 

where y res is the residual on subtracting the model from the data. Since in our cases the signal-to-noise is 
small in the channel data, this approach makes the estimate for measurement error robust and independent 
of errors in the model fit and any low-order residuals that are not represented in the model. 

The second model corresponds to the alternative hypothesis . This model contains the blackbody 
CMB term, the Maximally Smooth polynomial representing the foreground and, in addition, a template of 
the recombination lines that are expected to be present in the spectra for this hypothesis. As before, we 
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Fig. 14.— Bayes factors computed for data sets ‘a’ and ‘b’ that are mock observations generated with and 
without the recombination line component respectively. logio(BF) is plotted versus observing time assuming 
cryogenic cooled state-of-the-art 1 K receiver noise temperatures and a 128-element array of precision total- 
power spectrometers. Filled circles represent the median values of the FF’s computed separately for the 
two data sets. The horizontally-striped region shaded in red represents the range in Bayes factors obtained 
for the 100 mock observations in data set ‘a’ and the vertically-striped region shaded in green represents the 
range of Bayes factors obtained for data set ‘b\ 


proceed by fitting the data with the model, subtracting the best fit model from the data and computing 
the likelihood from the residual. We again make the reasonable assumption that if the data is completely 
represented by the model then the residual should be Gaussian random noise corresponding to measurement 
noise. The likelihood function P(D a \H 2 ) is given by: 


N 

P{Da\H 2 ) = n 

i=l 


_ Vres 2 [^] 
2°"2 




TTCTn 


( 21 ) 


where y re s 2 [i] is the residual on subtracting the alternative hypothesis model from the data, and variance o\ 
is half the Allan variance of this residual. Since the Allan variance gives the measurement noise independent 
of the model used to fit to the data, we expect cr^ and a\ to be the same. 


Having computed the likelihoods for both the models we arrive at the Bayes factor by computing the 
likelihood ratio: 


= P(Dg\H 2 ) 
P(D a \HoY 


( 22 ) 


We repeat this exercise with a separate data set c b’, which is a set of 100 spectra corresponding to mock 
observations with the same measurement noise variance as in data set ‘a’, but with no recombination lines 
in the synthetic spectra. We once again fit these spectra with the two models as before. We compute the 
likelihood functions for the null and alternative hypotheses given by P^D^Hq) and P(Db\H 2 ) respectively, 
under the reasonable assumption that when a model that completely represents the data, except for random 
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Fig. 15.— The confidence with which Bayes factor above unity might signify a detection, and simultaneously 
reject the detection to be a false positive, versus integration time. We assume cryogenic-cooled state-of-the- 
art 1 K receiver noise temperatures and a 128-element array of precision total-power spectrometers. 


noise, is subtracted from the data the residual must be the measurement noise. Thus the Bayes factor for 
the 100 mock observations that are generated assuming that recombination line component is absent is given 
by: 


BF b 


P(Db\H 2 ) 

P(D b \HoY 


(23) 


Fig. 14 shows a plot of logio {BF a ) and logio {BF^) versus integration time. With increasing integration 
time the measurement noise reduces and the two Bayes factors diverge, as expected. 


For any data set, the Bayes factor indicates the relative preference of the data for the two hypotheses. 
The median BF a rises above unity with increasing integration time, increasingly preferring the hypothesis 
H 2 over Hq. In contrast the median BF & drops below unity with increasing integration time, increasingly 
preferring the hypothesis that the recombination lines are absent in data set ‘b\ We may set a threshold at 
unity and examine the confidence with which this threshold might discriminate between observations that 
contain a recombination ripple to those in which this additive component is absent. At any integration time, 
the distribution of Bayes factors indicates the probability that BF a is above unity, which gives the confidence 
in a detection. For the same integration time, the distribution of BF & yields the probability of obtaining a 
value above unity, which gives the likelihood of a false positive. Using the multiple data sets ‘a’ and ‘b’ we 
have computed the run of the confidence in detection (and simultaneous rejection of false positives) versus 
integration time: this is depicted in Fig. p~5[ 

If we adopt the Bayes factor as a detection statistic, detection of the recombination ripple with 68% 
confidence using uncooled receivers requires observing with about 52 x 10 8 antenna secs, or equivalently about 
60,000 antenna days. An array of 128 singly-polarized total-power spectrometer elements with uncooled 
receivers would require about 1.3 years for detection with this confidence. Using cryogenically cooled state- 
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of-the art receivers, the observing time for such detection with 68% confidence is 2.56 x 10 8 antenna secs or, 
equivalently, 23 observing days with an array of 128 spectral radiometers deploying such cooled receivers. 
90% confidence detection requires observing with 255 days with such an array, and 95% confidence detection 
requires about 440 days. 

Fig. [16] shows the confidence in the detection of cosmological recombination lines in synthetic sky spectra 
as a function of integration time, estimated using both the Bayes factor method as well as using the fit to 
the amplitude of the recombination ripple (Section |5.1[ ). The Bayes factor approach to detection is simplest 
in that it only asks whether in an observed spectrum the lines as predicted by theory are present or totally 
absent. The method that derives an amplitude for the ripple, as estimated by a best fit to a scaling factor 
for a template of the expected ripple, attempts to ask a more informative question in that it attempts to 
evaluate the amplitude of the recombination line ripple. Unsurprisingly, the Bayes factor method has greater 
confidence in its answer for any integration time. 



Fig. 16.— A comparison of the confidence in detecting cosmological recombination lines in sky spectra 
using Bayes Factors and that given by a comparison of the best-fit ripple amplitude to a threshold. The 
confidence values are shown versus integration time assuming cryogenic cooled state-of-the-art 1 K receiver 
noise temperatures and a 128-element array of precision total-power spectrometers. 


6. Summary 

We demonstrated that it is in principle feasible with current day technology to experimentally detect 
the cosmological Hydrogen and Helium recombination lines at low frequencies, although the lines are embed¬ 
ded in a foreground that is about nine orders of magnitude brighter with a priori unknown precise spectral 
shape. The recombination radiation has a smooth component that is difficult to distinguish from the smooth 
foregrounds; however, the recombination radiation also has a unique ripply component that may be distin¬ 
guished from the foregrounds and instrumental effect. We estimate the amplitude of the ripply signal, the 
instrument noise arising from receivers and foregrounds, and estimate the signal-to-noise in detections using 
ground based spectrometers. Detection may ideally be attempted using an octave band in the 2-6 GHz 
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window; an octave band would have a spectral segment of the recombination ripple with sufficient struc¬ 
tural complexity so as to be distinguishable from the relatively smoother foreground. We have developed an 
algorithm to detect the recombination line ripple by foreground modeling that models the foreground as a 
Maximally Smooth polynomial, which we define enforcing no zero crossings in derivatives of order n > 2. A 
similar approach may be applicable to modeling and detecting the global 21-cm signal at lower frequencies. 

We then evaluate the confidence with which these cosmological recombination lines may be isolated 
using the aforementioned algorithm and estimate the integration times required for detection with varying 
degrees of confidence. In its simplest form, a detection is testing an observed spectral segment addressing 
the question of whether or not the theoretically motivated spectral ripple is present in the sky spectrum or 
absent. To answer this question with 90% confidence requires an integration time of 32,640 antenna-days 
with a total-power spectral radiometer with cryogenically cooled state-of-the art receivers; however, with 
today’s best uncooled receivers the integration time for such a detection increases to as much as 660,000 
antenna days. This makes a compelling case for the use of cryogenically cooled receivers. Moreover, since the 
antenna element for such all-sky global signals does not gain from antenna directivity, very small antenna 
elements of centimeter dimensions may be deployed, perhaps as a compact cluster of antenna elements housed 
in a moderate-size dewar. Assuming a 128 element array followed by independent spectral radiometers, a 
90% confidence detection may be achieved in 255 days integration time. 

APSERs^j— Array of Precision Spectrometers for the Epoch of RecombinAtion—is an experimental 
venture at the Raman Research Institut^J India, with the science goal of detecting spectral ripples from 
the epochs of Hydrogen and Helium recombination. On completion, it is to be an array of 128 miniature 
radio telescopes operating in the octave band of 3-6 GHz with antenna elements and receivers that are 
custom designed for all-sky spectral measurements with the required sensitivity to experimentally detect 
these cosmological recombination lines. APSERa will be deployed at a radio quiet site closer to the geographic 
poles to avoid radio frequency interference in these bands from downlinks of geostationary satellites. 


7. Future Investigations 

Ripple like spectral signatures arising from the epochs of Hydrogen and Helium recombination are 
predicted to appear as additive distortions to the sky spectrum at a level ^9 orders of magnitude weaker 
than the total brightness of the radio sky. We have demonstrated that it is in principle possible to detect 
these signals. Due to their small amplitude, any small departure in the spectral radiometer from the ideal 
behavior assumed herein could critically affect the detection likelihood. Studies of the effects of potential 
instrumental non-idealities, and confusion arising from additive astrophysical contaminants, which have not 
been considered here, will form the next step in our future work. The key investigations are listed below. 

• The spectral radiometer radio telescope used to detect the weak cosmological recombination lines 
must be capable of achieving the required spurious-free dynamic range so as to recover the embedded 
cosmological signal in the measurement. This implies that the receiver configuration and calibration 
techniques must be specifically designed to minimize systematics and that any residual systematics are 
either well characterized or emendable to accurate modeling so as to enable unambiguous detection of 


6 http://www.rri.res.in/apsera/ 
7 http://www.rri.res.in/ 



the weak signal. A future work is the design of a suitable receiver configuration, observing strategies, 
gain stability and calibration methods to achieve the required accuracy in the measurement. 

In our simulations, we have generated mock sky spectra as observed by an ideal frequency independent 
antenna beam. The spectral response of such an antenna is maximally smooth if all sources in the 
beam have spectra that are complete monotones. However, designing and fabricating an antenna that 
is frequency independent over an octave bandwidth is non-trivial. Frequency dependence in the beam 
pattern and its sidelobes may result in a response that is no longer smooth and may confuse the 
detection of cosmological recombination. A study of this mode coupling of sky structure into spectral 
structure in a spectral radiometer, which depends on the type of frequency dependence in the beam 
pattern, would lead to design tolerances on the antenna element for this detection experiment. 

The all-sky ‘ripply’ signal arising from the epochs of Hydrogen and Helium recombination is inher¬ 
ently unpolarized and this spectral feature may be detected with an antenna that responds to any 
single polarization mode, either linear or circular. However, the foreground synchrotron emission from 
extragalactic sources as well as Galactic emission is linearly polarized and Faraday rotation during 
the line-of-sight propagation results in the received polarization position angle varying with observing 
frequency. This causes linearly polarized sources to appear with a ‘ripply’ spectral structure in the 
response of linearly polarized antennas. One possible way of eliminating beam asymmetries and the 
effect of beam rotation across the sky is to rotate the array itself. A study of this mode coupling of 
polarized sky emission to spectral structure is a future study that will lead to design tolerances on the 
polarization properties of the antenna element and on the polarization calibration. The potential of 
using the unpolarised nature of the cosmological signal to discriminate it from polarised foregrounds 
is another aspect to be explored. 


The spectral template of the additive ripple-like feature from cosmological Hydrogen and Helium 
recombination has a fairly accurate theoretical prediction. The near quasi-sinusoidal nature of this 
signal acts as a fingerprint providing a means to distinguish it from other weak cosmological signals. 


Two such cosmological signals are the (i and y distortions of the CMB (Zeldovich & Sunyaev 1969 


Sunyaev & Zeldovich 1970). Even within the standard cosmological model, these distortions are created 


at amplitudes that are typically larger than the cosmological recombination radiation (e.g., Chluba & 


Sunyaev 2012 Sunyaev & Khatri 2013 Chluba 2013 Andre et al. 2014). A large average Compton 
y distortions, at the level of y ~ 10~ 7 — 10 —6 , is expected from reionization and structure formation 
( Sunyaev fc Zeldovich|1972 Hu et al.|1994| Cen fc Ostriker|1999 Oh et al.|2003 ) and unresolved clusters 
and filaments (Refregier & et al. 2000 Miniati et al. 2000; Zhang et al. 2004). At low frequencies, this 
type of distortion only leads to a weak frequency dependent signal, T(v) ~ —2^/(1 — x 2 /12) with 
x « 0.018 [z//GHz], making it is less problematic. For a /i distortion, the low frequency spectrum 
shows much richer structure, especially when varying time-dependence of the energy release mechanism 
(Fig. 14, Chluba & Sunyaev 2012). Each of these introduce spectral features in the CMB with the 
/i distortion peaking at about 1 GHz. In an experiment optimized to detect ripples originating from 
the epochs of Hydrogen and Helium recombination it would be interesting to study whether /i and y 
distortions of the CMB are a possible source of confusion or would themselves be detected as a positive 


by-product. Similarly, the anomalous microwave emission from spinning dust (see Draine & Lazarian 


1998; Ali-Hai moud et al.|2009 Planck Collaboration et al.|2014 ) will add another spectral dependences 


to the low-frequency regime, which at this point we have not investigated. 
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A. The Completely Monotone (CM) function 


Mathematically, a function f{pc) is said to be a complete monotone (CM) if for all values of x in the 
interval 0 < x < oo, (—l) n x d > 0 for every integer n > 0. The method by which a CM function of any 
said order may be constructed is best explained by an example, such as that for a cubic function as given 
below. We construct a cubic function f(pc) which is CM in the range defined by [xo,x m \: 


f(x) = Po - Pi(x - X 0 ) + p 2 (x - x 0 ) 2 - p 3 (x - x 0 ) 3 , (Al) 

where Pi for all 0 < i < 3 are the coefficients that are constrained so that (—1)” x d d ^„ > 0 for every 
integer n > 0 and for all x in [xq . x m ] . The third derivative of f(x) is 


d 3 /Qr) 

da: 3 


-6p 3 - 


(A2) 


To constrain f(x) to be CM we require that (—l) 3 x d > 0, or that p 3 > 0. Let p 3 = <23, where <23 is a 
positive value. The second derivative of f{pc) is 

d 2 f(x) , N 

dx 2 = 2p 2 - 6ps{x - Xo) 

= 2p 2 - 6a 3 (x - x 0 ) 

Once again to constrain f(pc) to be CM we require that (—l) 2 x d > 0, which leads to the condition that 
2p 2 — 6a 3 (x — xo) > 0. This is satisfied if P 2 is selected such that 2p2 — 6a 3 (x rn — xo) > 0, which ensures that 
the criterion is satisfied for all x in [xo,x m \. Therefore, we express P 2 as P 2 = <22 + 3a3 (x max — xo), where 
a 2 > 0. The first derivative of f(x) is 


df(x) 

dx 


-Pi + 2p 2 (x - xo) - 3 ps(x - xo) 2 

-Pi + 2{a 2 + 3 a 3 (x m - x 0 )}(x - x 0 ) - 3 a 3 (x - x 0 ) 2 


To constrain f(x) to be CM, we once again require that (—1) x d ^^ > 0. This requires that we determine 
Pi such that — pi + 2 (a 2 + 6as{x m — xo)}(x m — xo) + 3(a3 )(x m — xo) 2 > 0, ensuring that the criteria is 
satisfied for all x in [xo,x m \. To guarantee this we express pi as pi = ai + 2a2(x rn — xo) + ?>as{x rn — xo) 2 , 
where a\ > 0. Since the constant term po represents a vertical translation of the function and does not affect 
the shape—the smoothness of the function itself—we allow it to be a free parameter < 20 . In summary, we 
describe the cubic CM function to have the form 


f(x) =a 0 - {ai + 2 a 2 (x m -x 0 )+ 3 a 3 (x m - x 0 ) 2 } (x - a; 0 ) + {a 2 + 3a 3 (x m - x 0 )j(x - x 0 ) 2 - a 3 (x - x 0 ) 3 , (A3) 

where all of the coefficients ao, ai, a 2 and a 3 are non-negative. In this functional form, setting a 3 = 0 reduces 
the function to CM of order 2 and if we further set <22 = 0 we have a CM polynomial of order unity. This 
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formulation allows a successive approximation to data by first fitting for a CM polynomial of order unity as 
a Taylor expansion at the reference point x = xq to optimize for coefficients ao and a\ , then adding a higher 
order term with initial guess — 0 and optimizing all three coefficients, and finally adding a cubic term 
with initial guess as = 0 and optimizing all four coefficients. All coefficients are constrained to be positive 
by expressing them in the form ai = 10 6i , where the are allowed to be real numbers. 

Thus a CM polynomial of arbitrary order n may be written in the form 

n n—i 

f(x) = a 0 + yjfe l )'(.!• - zo)'{y ^aj+jCj +J (x m - x 0 ) J }, (A4) 

i=1 3 =0 

where denotes the binomial coefficient n\/{k\(n — fc)!} and the constant term ao is left unconstrained so 
as to allow arbitrary vertical translations. 
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